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Abstract. The very high energy Galactic 7 -ray sky is partially opaque in the (0.1 — 10) PeV 
energy range. In the light of the recently detected high energy neutrino flux by IceGube, a 
comparable very high energy 7 -ray flux is expected in any scenario with a sizable Galactic 
contribution to the neutrino flux. Here we elaborate on the peculiar energy and anisotropy 
features imposed upon these very high energy 7 -rays by the absorption on the cosmic mi¬ 
crowave background photons and Galactic interstellar light. As a notable application of our 
considerations, we study the prospects of probing the PeV-scale decaying DM scenario, pro¬ 
posed as a possible source of IceGube neutrinos, by extensive air shower (EAS) cosmic ray 
experiments. In particular, we show that anisotropy measurements at EAS experiments are 
already sensitive to tdm ~ 0 ( 10 ^^) s and future measurements, using better gamma/hadron 
separation, can improve the limit significantly. 
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1 Introduction 

The very high energy (VHE) 7 -ray band, loosely dehned as the one above 0(100) GeV, is 
peculiar in many respects. For instance, on the detection side, it is at the upper edge of the 
energy range which can be probed from space via direct techniques, and can be effectively 
studied only by ground-based telescopes. An important theoretical feature in the study of 
VHE 7 -rays is that absorption is relevant, i.e. the extragalactic sky is not transparent to VHE 
7 -rays due to the pair-production absorption onto the Extragalactic Background Light (EBL) 
and the cosmic microwave background (CMB). Actually, for photon energies of 0(100) TeV 
even the Galactic diffuse radiation held starts playing a role, while at ~ PeV energies the 
mean-free path of 7 -rays due to the absorption on GMB photon bath approaches 0(10) kpc, 
comparable to the distance to Galactic Genter (GG). Hence, for VHE 7 -rays even the Galactic 
sky becomes partially opaque. While the opacity of the extragalactic sky is routinely taken 
into account in the analysis of cosmologically distant active galactic nuclei, both by the Fermi- 
LAT space telescope [1] and by the ground based instruments, see e.g. [2], the absorption at 
Galactic scale is usually neglected, although the theoretical importance of the phenomenon 
has been acknowledged in the past, see e.g. [3]. 

The most obvious reason for the relative lack of interest in this phenomenon is that 
the “PeVatron” window in 7 -ray astrophysics has yet to be opened, with current telescopes 
running out of statistics below 0(100) TeV (even for the most powerful sources) to be sensitive 
to spectral distortions. Another reason may simply be that PeVatrons (in gamma-rays) in our 
Galaxy are expected to be rare, if present at all, requiring extreme acceleration parameters. 
The situation has however changed in the last couple of years, following the discovery by the 
IceGube collaboration of a diffuse neutrino flux in the 30 TeV-2 PeV range [4-6]. Virtually in 
any conceivable model for its origin, one expects an associated 7 -ray flux with similar energy 
budget in the ~ PeV energy range. The main unknown concerns the distance at which the 
sources of these events are located: if the neutrino flux is due to a population at cosmological 
distances, the absorption should be so severe that the initial VHE flux would cascade down 
below 0(100) GeV, contributing to the residual isotropic background precisely measured by 
Fermi [7]. However, if it is due to Galactic or nearby extragalactic sources (see e.g. [8-11]), 
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the associated gamma-ray flux should still be detectable at VHE, albeit with a signiflcantly 
suppressed spectrum. This has been noted soon after the discovery, see for instance [12] 
or [13], and current 7 -ray constraints are one argument disfavoring close-by discrete Galactic 
sources for the majority of the events. 

The calculation of both the expected signal and the observational constraints is however 
more involved in the case of truly diffuse local sources associated with the IceCube data, such 
as in an astrophysical origin in the Galactic halo [14], or a decaying Dark Matter (DDM) 
origin, see [15-19] (see also [20-23] for some DM related interpretations of the IceGube data). 
On the one hand, the expected signal is usually anisotropic, at very least due to our off-set 
position with respect to the center of the Galactic Halo, so its detailed calculation requires 
at least a 2D modeling of the problem, and multiple numerical integrations, a complication 
that typically is not fully taken into account even in the most recent calculations [24]. On 
the other hand the observational bounds (the most constraining ones being the one published 
by GASA-MIA [25] and those reported in a proceeding by KASGADE [26]) are derived: a) 
based on observations of a limited portion of the sky; b) typically assuming a isotropic 7 -ray 
flux. 

The hypothesis b) is clearly incorrect, since even an incoming isotropic flux would ac¬ 
quire an anisotropy due to the anisotropic absorption. The hypothesis a) means that a proper 
use of these constraints would require to re-run ad hoc analyses by the collaborations, based 
on specific energy and angular-dependent templates, to be convoluted with the detector char¬ 
acteristics. This is a pity, since as we have argued in [16], for an important class of scenarios 
like the DDM ones, this kind of data are what currently comes the closest to an independent 
test of the hypothesis 

Although the important role of EAS probes of this scenario has been discussed in the 
past (see e.g. [12, 16, 24]), here we revisit the calculation of the expected 7 -ray flux, with a 
triple goal: i) To estimate more precisely the spectral and angular shape of a DDM signal, 
with state of the art treatment for the primary 7 -ray absorption and the inverse Gompton 
component, ii) To point out that due to the generically anisotropic nature of the VHE 7 -ray 
component, even detectors with little or without gamma-hadron rejection capability should 
be able to put constraints on these contributions based merely on the expected anisotropy. 
Hi) To motivate experimental collaborations to specifically constrain some angular-energy 
templates, to optimize their constraining power for specific models. For the DDM case, for 
instance, an intermediate step in this direction would be to derive (energy-dependent) bounds 
in coronas around the GG, in Galactic cylindrical coordinates. We also discuss the greatly 
improved potential of detectors with significant hadron vs. 7 -ray rejection capabilities. 

As a case study we consider throughout this paper the peculiarities of the expected 
7 -ray flux from DDM. Yet, similar considerations would apply to any other Galactic diffuse 
flux model (a few examples have been listed e.g. in [27]). We will consider both the prompt 
7 -ray flux and the flux from inverse Gompton (IG) scattering of off the ambient photon 
bath. Both contributions would be present also in other diffuse flux models: in the commonly 
considered astrophysical hadronic production of neutrinos, they are associated to prompt 7 - 

^Note that a superficial look at Fermi-LAT isotropic gamma-ray background (IGRB) results might suggest 
that they are already very constraining, notably thanks to the few high-energy points. We stress here that they 
are unfortunately not robust with respect to the IGRB extraction procedure from the Extragalactic gamma- 
ray background (EGB). In more technical terms, the determination of the IGRB by the Fermi-LAT team does 
not take into account the uncertainty in the subtraction of the point sources contribution (very uncertain at 
high energies, relying on an extrapolation), which would constitute the dominant source of uncertainty in the 
last IGRB points. If the EGB is conservatively used, Instead, the bounds are degraded, as can be seen in [24]. 
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rays from vr^ decay and from vr^ decay, respectively. The specific features studied in this 
paper, mainly the inherent anisotropy due to absorption at Galactic scale and the peculiar 
profile of IC flux due to diffusion/losses of PeV , would similarly provide powerful diagnostic 
tools in probing alternative diffuse flux models. The only differences would be in the initial 
spectra and the geometric distribution of the source term in the Galaxy/Galactic halo. 

This article is structured as follows: in section 2 we describe the peculiar energy-angular 
dependence of the 7 -ray flux absorption, and our computation of the 7 -ray opacity. In sec¬ 
tion 3 we compare the expected 7 -ray flux from DDM with current constraints from EAS 
experiments, as well as the diagnostic power of forthcoming experiments. The two compo¬ 
nents of the 7 -ray flux from DDM, prompt and IG flux, are discussed in sections 3.1 and 3.2, 
respectively. Section 4 is devoted to the discussion of the expected anisotropy of the total 
cosmic ray flux. Finally, in section 5 we conclude. 


2 Absorption of 7 -rays at Galactic scale 


The 7 -ray flux in the approximate range 10“^ -7 10^ PeV will suffer attenuation in the Galaxy 
due to the pair production 77 —)■ e^e"*" process onto photon baths: at the lower energies, 
starlight (SL) and infrared (IR) photons constitute important targets (mostly however for 
directions towards the inner Galaxy), while at ~ PeV energies and above the homogeneous 
cosmic microwave background (GMB) is dominant. In the following we calculate the optical 
depth for both GMB and SL+IR, for different incoming directions and energies. 

For the technically simpler case of pair production on GMB photons, the optical depth 
for photons of energy coming from a source at distance L can be calculated as (here and 
in the following, we use natural units with c = /cs = 1 ) 


= L 11 a^^{E^,e)ncMB{e) 


where is the pair production cross section given by 


1 — cos 6 


sin 0 d 0 de , 


( 2 . 1 ) 


^77 = ^^(1 - (3 - /3") In ( [^ ) - 2/3(2 - (5^) 


( 2 . 2 ) 


where 


/3 = 


, and 


s = 


eEr^ 

2ml 


(1 — cos 9 ), 


(2.3) 


a is the fine-structure constant, rUe the electron mass and 9 is the angle between the momenta 
of photons. The ncMB(£) is the differential number density of GMB photons given by 
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where Tcmb = 2.348 x 10““^ eV. By changing variable e —)• £c, where £c = ^j£E.y{l — cos0)/2 
is the photon center of momentum energy, and performing the integral on 0 , the expression 
for can be reduced to a single integral to be performed numerically, namely 
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Figure la shows the as function of Er^ for three different values of L = 4 kpc, 8.3 

kpc and 20 kpc. As can be seen, for a source of 7 -ray at Galactic center (GC), at about 
L = 8.3 kpc, the absorption is ~ 70% at ~ 2 PeV. Figure lb shows the contour plot of 
exp[—as function of photon energy E^ and source distance L. 



Figure 1. Plot of the absorption of 7 -rays on CMB photons. Panel (a): for a source at distance 
L = 4 kpc, 8.3 kpc and 20 kpc. Panel (b): 2D density plot of exp[—as function of L and E^. 


The optical depth due to pair production on the SL+IR photon bath can be calculated 
similarly to eq. ( 2 . 1 ), with the extra complication that the integral along the line of sight is 
non-trivial, since the photon bath number density nsL+lR also depends on position x, and 
the optical depth also depends on the Galactic coordinates {b,l). In the approximation that 
the photon field is inhomogeneous but isotropic one can write 

= j ds jj a^^{E^,e)nsL+iR[£,^{s,b,l )]-— ^^sinMede, (2.6) 

where the line-of-sight parameter s is related e.g. to the cylindrical coordinates (r,z), with 
the origin at the GG, by 


r = y cos^ 6 — 2 si ?0 cos 6 cos/ and z = ssinb, (2-7) 

where Rq ~ 8.3kpc is the distance of the Sun to the GG. The number densities of SL and 
IR photons have been extracted from the GALPROP code [28] and their energy densities for 
some representative positions are plotted in figure 2. Obviously, the GMB radiation field is 
homogenous and thus pervades the whole Galaxy uniformly, while the SL and IR components 
of radiation field are clearly position dependent: larger at GG and in the Galactic disk, 
decreasing rapidly by moving perpendicularly from Galactic disk, along the 2 ; direction. 

The optical depth due to SL+IR photon bath for two different distances and various 
directions are shown in figure 3. It is clear that the absorption effect is relevant around 
energies of 0(100) TeV, but only for directions towards the inner Galaxy (6 ~ ~ 0). The 

calculated optical depths in this section are consistent with the results reported in [3]. The 
effect of the total opacity of Galactic medium (i.e., will be discussed 

in the following sections. 


4 












Figure 2. The energy density of ISRF (including starlight, IR and CMB), extracted from GAL- 
PROP [28], at three different positions in our Galaxy: the dotted curves are for (r, z) = (0, 0), that is 
GG; the solid curves are for (r, z) = (8.3,0) kpc, that is the Sun position, and the dot-dashed curves 
are for (r, z) = (8.3, 5) kpc. 



Figure 3. Plot of the absorption of y-rays on SL+IR photons, for a source at distance L = 8.3 kpc 
(solid curves) and 20 kpc (dashed curves) for various directions. 
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3 Expected fluxes from DDM and constraints from EAS experiments 


3.1 Prompt component 

The prompt component of the Galactic 7 -ray flux from DM decay in the direction ( 6 , 1) is 
given by 


dE^ 


{Ej, 6, 1) 


- - ^r ds , 

dvr ttidm t dm dE^ Jq 


(3.1) 


where moM and tdm are respectively the DM mass and lifetime, and r^yy is the total optical 
depth, ph is the density prohle of DM particles in our Galaxy as a function of radial distance 
(in spherical coordinates) from the Galactic center, g. For our hducial model we adopt a 
Navarro-Frenk-White density prohle [29] 


Phig) ^ 


Ph 

g/rdl + q/tcY ’ 


(3.2) 


where Tc — 24kpc is the critical radius and p^ = 0.18 GeVcm“^, which yields a DM density 
at the Solar System pq = 0.39 GeV cm“^ [30]. The line-of-sight integration parameter s is 
related to radial distance g via 


g{s, b,l) = \J+ Rq — 2sRq cos b cos I. (3.3) 

The dNy /dEy is the energy spectrum of photons produced in the decay of a DM particle (here 
obtained from the PYTHIA 8.2 [31], including the weak gauge boson radiation corrections as 
from [32]). To illustrate the typical spectra from DM decay, in hgure 4 we plot the EydN/dEy 
for various decay channels of a DM particle with moM = 4 PeV. In a specihc model of the DM 
(i.e, specihc decay channels with branching ratios determined by the model) the spectrum of 
7 -ray can be obtained by the appropriate weighting of the spectra in hgure 4. In this paper 
we adopt the scenario introduced in [33] where the heavy DM particle is a sterile neutrino 
with mass ~ 4 PeV and lifetime ~ 10^® s, with the branching ratios of the decay channels 
given by 

Bv{e^W^) = 2Br(yZ) = 2Br{^Y!h) = \Ueif , (3.4) 

where U^i are the the elements of the neutrino mixing matrix (for details see [33]). It is shown 
in [16] that this scenario provides a reasonable ht to the energy distribution of IceGube 
neutrino data. However, let us emphasize that these choices of branching ratios and decay 
channels are not extremely constrained. In fact, as discussed in [15], any model with a sizable 
branching ratio into hard (leptonic) channels, with the remaining (even dominant) branching 
ratio into soft (hadronic and gauge bosons) channels, would provide decent ht to the IceGube 
data. 


3.2 Including the inverse Compton component 

An additional component of the 7 -ray hux comes from the inverse Gompton (IG) scattering 
of the electrons and positrons from DM decay, up-scattering mostly the GMB photons, which 
writes 
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Figure 4. Spectrum of the 7 -ray yield in various decay channels of a DM particle with todm = 4 PeV. 
The red curve for DM —>■ quarks, shows the average spectra for the DM decay to all the quark flavors. 
The spectra are obtained via PYTHIA 8.2 [31]. 


where g is given in eq. (3.3), Pic is the IC power due to up-scattering on different photon 
backgrounds and dUe/dEe is the differential number density of e"*" plus e~ at steady state. 
Although the IC flux reported in this article is calculated by taking into account the spatial- 
dependent nature of energy losses and the effect of spatial diffusion (see appendix A for the 
details of the calculation), in order to grasp the main features of IC flux, in the following we 
pursue a simplified version of the calculation. At the energies of interest here the transport of 
electrons and positrons in our Galaxy is determined almost exclusively by the energy losses. 
Also, one realizes that for directions close to the Galactic plane and for realistic values of 
the Galactic magnetic field (i.e., ~ few /rG, with a profile that increases towards the inner 
Galaxy) synchrotron emission is the dominant energy loss mechanism, simply because the 
synchrotron emission is always quadratic in the electron energy and does not suffer the Klein- 
Nishina suppression of IC on SL and IR photon baths. Also, at high energies the IC power 
Pic is almost exclusively due to up-scattering of the CMB photons, and thus independent of 
the position. The position dependence of the energy loss coefficient, b = —dEe/dt, is more 
involved and traces the Galactic magnetic field profile. However, in the approximation in 
which the thin gaseous disk of the Galaxy is embedded in a thick diffusive halo permeated 
by a constant magnetic field, the loss coefficient b is independent of the position. In this 
approximation (which we checked to be accurate whenever the IC signal is non-negligible, see 
appendix A for details), one can write 


dne 
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(3.6) 
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and the total 7 -ray flux (i.e., sum of the prompt and IC components) can hence be written 
as 
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where 


dA^., 
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ANe/AEf, is the energy spectrum from DM decay, obtained via PYTHIA 8.2 [31]. Pic can 
be calculated straightforwardly as reported in appendix A. Yet, the energy loss coefficient b 
still depends on the poorly known value of the magnetic held, Shalo permeating the thick 
halo which extends to several kpc away from the disk, and for the lack of better information 
we approximated it as constant. 

In figure 5 we show the 7 -ray flux from DM decay, assuming moM = 4 PeV and tdm = 
10 ^® s (chosen to be close to best-fit parameters; the flux scales inversely with ut-dm and tdm) 
and for the decay pattern with the branching ratios of the decay channels given by eq. (3.4), for 
different directions in Galactic coordinates. The solid curves depict the prompt flux, eq. (3.1), 
from GC (red, top), anti-GG (blue, bottom) and Galactic Pole (orange, intermediate). In each 
of these curves the dot-dashed curve deviating from the solid curve at higher energies shows 
the flux neglecting the absorption of 7 -rays discussed in section 2. When comparing the 
expected 7 -ray flux from DDM with the experimental bounds, the importance of accounting 
properly for the absorption of 7 -ray on GMB photons is manifest, particularly at high energy. 
The dashed curves in figure 5 show the IG flux: the red (blue) dashed curves are the IG flux 
form GG (anti-GG) direction. The orange dashed curve shows the IG flux from the Galactic 
pole direction, with the assumption that the Galactic magnetic held only consists of the (thin 
disk) regular held, given in eq. (A. 5); i.e., i?halo = 0- The cyan, black and green dashed curves 
show the IG flux from the Galactic pole within the assumptions i?haio = 0.5, 1 and 2 /xG, 
respectively. Finally, the green and brown bar lines with arrows show, respectively, the upper 
limits on the 7 -ray flux inferred by GASA-MIA [25] and KASGADE [26] experiments. 

Let ns elaborate on the various IG fluxes shown in figure 5: the IG component is clearly 
sub-leading with respect to the prompt emission for directions along the Galactic plane (note 
the red and blue dashed curves). However, this is not necessarily the case for the IG flux 
from the Galactic poles. The enhancement of the IG flux from the Galactic pole direction 
originates from the fact that for vertical directions the b coefficient drops fasters than the 
DM density along the line of sight integration of eq. (3.7). This enhancement is sizable 
if one assumes that the magnetic held exponentially decreases for vertical directions (with 
the profile of eq. (A. 5))— dashed orange curve in figure 5— so that the IG flux can become 
comparable to the prompt flux towards the Galactic pole. However, as we have mentioned 
earlier, it is realistic to assume that a non-zero magnetic held permeates a thick halo to 
large distances, as consistent with the assumption that a charged cosmic ray population still 
propagate diffusively in a region several kpc away from the disk. The constant Hhaio is a 
toy-model representation of this held, and its effect on the IG flux can be seen by the cyan, 
black and green dashed curves. In all cases, the emission is suppressed with respect to the 
“unmagnetized halo” situation. The reason is that a growing Hhaio leads to a larger energy 
loss coefficient b, and thus more suppressed IG flux, since a growing fraction of energy is 
channeled into synchrotron. In conclusion, the IG flux from directions close to the Galactic 
plane (low- 6 ) is quite robustly predicted to be small. The exact value of IG flux towards the 
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Figure 5. The 7 -ray flux from DM decay from various directions, with todm = 4 PeV and tdm = 
10^® s, and branching ratios of decay channels given by eq. (3.4). The solid curves are shows the prompt 
flux, including the absorption of 7 -rays, while in the dot-dashed curves the absorption is neglected. 
The dashed curves show the IC flux, for various assumptions for the constant halo magnetic held, 
43haio, possibly pervading the thick diffusive halo of the Galaxy up to large distances. The green 
and brown bar lines show the upper bound on 7 -ray flux from CASA-MIA [25] and KASCADE [26] , 
respectively. 


Galactic poles is hard to predict due to the uncertain thickness and i?-field strength of the 
magnetized halo, with the orange dashed curves in figure 5 providing a reasonable upper limit 
to this uncertain component. 

It is worth noting that the CASA-MIA and KASCADE experiments would have al¬ 
ready probed interesting parameter space for DM models, if they had accumulated significant 
exposure towards inner Galaxy, e.g. if they had been located in the Southern hemisphere. 
Unfortunately, their acceptance mostly peaks in regions far away from the GC and hence 
they would have been exposed to more modest fluxes, comparable to the orange curve in 
figure 5, insufficient to test the model even for optimistic IC expectations. To illustrate this 
point, in the following we briefly describe some notions on the geometrical acceptance of 
EAS experiments. An EAS is often classified as 7 -like event, as opposed to a hadronic-like 
event, based on a significantly poorer muon content of the former shower with respect to the 
latter (at a fixed primary energy). Only for events which are not too inclined with respect 
to the vertical this separation can be done meaningfully, thus imposing a cut on maximum 
zenith angle of the shower. Assuming that the detector is continuously operational (i.e., the 
acceptance is uniform with respect to azimuth, or right ascension in equatorial coordinate), 
the geometrical acceptance efficiency u) of an EAS experiment located at the latitude A as 
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function of declination S, can be written as [34] 


uj(S) oc cos A cos S sin Om + Om sin A sin S , 


(3.9) 


where 




0 c > 1 

vr C < “1 

arccos(C) —1 < C < 1 


cos — sin A sin 6 

and C =-r- 7 - 

cos A cos 0 


where 0m is the maximum zenith angle acceptance of the detector. For the CASA-MIA ex¬ 
periment (A, 0m) = (40.2°, 60°) [25], and for KASCADE (A, 0m) = (49°, 35°) [26]. The blue 
shaded regions in hgures 6 a and 6 b show the estimated acceptance of CASA-MIA and KAS¬ 
CADE experiments, respectively, in Galactic coordinates. Darker blue regions correspond to 
higher effective exposure to a 7 -ray flux. The superimposed red curves on the plots are the 
contours of at Ej = 1 PeV, with values 2 x 10“^^ down to 3 x 10“^^, respectively 

from the inner to outer circles (in units TeV cm“^ s“^ sr“^), from the decay of DM with 
^DM = 4 PeV, tdm = 10^^ s and branching ratios given in eq. (3.4). From hgure 5, it follows 
that the upper limits on E^d^y/dEj at Ej = 1 PeV, are ~ 2 x 10“^^ and ~ 4 x 10“^^ respec¬ 
tively from CASA-MIA and KASCADE experiments. However, these upper limits apply to 
the dark blue regions of hgure 6 and, as can be seen, the limits relax by moving to the regions 
where the 7 -ray hux from DM increases. In fact in the regions where these experiments are 
mostly sensitive, the expected hux from decaying DM is ~ 3 X 10 TeV cm ^ s ^ sr 
which is almost one order of magnitude below the KASCADE upper limit. 


4 Anisotropy 


Despite the fact that current EAS bounds are not yet constraining enough for the DDM 
explanations of IceCube events, the interesting parameter space appears within reach. Even 
relatively modest optimizations of current sensitivities might thus prove crucial. In fact, the 
main reason for the degradation of the bounds discussed in the previous section relates to 
the incorrect assumption that the gamma-ray hux is isotropic. In this section, we discuss 
to what extent one may turn that weakness into an opportunity, suggesting that anisotropy 
studies alone, even without shower property discrimination capabilities, might contribute to 
the constraints. EAS experiments in fact routinely measure cosmic ray anisotropy, albeit 
often only in terms of some “partial estimator” like the dipolar anisotropy (averaged with 
respect to right ascension). Let us dehne a characteristic “gamma-ray induced anisotropy” as 


a.y = 


dE~/ dE^ 

dE 


anti—GC 


(4.1) 


where d^cR/dFl is the total cosmic ray hux, taken from [35]. The anisotropy variable as 
dehned in eq. (4.1) mainly arises from the prompt hux and the contribution of IC hux is 
negligible, not only because the IC is sub-leading but also since it is expected to be relatively 
more isotropic. An immediate constraint on DM lifetime can be obtained by requiring that 
a-y does not exceed the observed total anisotropy in cosmic rays, a. In practice, by requiring 
that in no energy bin exceeds by more than two sigma the measured value of a we can 
obtain a conservative bound on the DM lifetime as tdm > 2.5 x 10^^ s. The power of this 
observable is due to the fact that the intrinsic anisotropy in charged cosmic rays is at the 
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(a) CASA-MIA 



Figure 6. The Mollweide plot of the efficiency of EAS experiments (a) CASA-MIA and (b) 
KASCADE, in Galactic coordinates, discussed in eq. (3.9). The red curves show the contours of 
/dEry Sut Ery = 1 PcV, wfth valucs 2 X 10“^^ down to 3 x 10“^^, respectively from the inner to 
outer circles (in the unit TeV cm“^ s“^ sr“^), from the decay of DM with rriDM = 4 PeV, tdm = 10^® s 
and branching ratios given in eq. (3.4). 


level of 10“^ -b 10“^, while a much larger (by two to three orders of magnitude!) relative 
anisotropy in gamma-rays is expected, at very least due to the off-center position of the Sun 
in the Galaxy. This means that, despite the fact that gamma-rays only constitute a small 
fraction of the overall CR flux at 0.1 — 1 PeV energies, in the anisotropy observable one can 
beneht from a larger signal to noise ratio. Accounting for absorption, however, suppresses 
the gamma-ray anisotropy, since pair-production is more severe in the GC direction than the 
anti-GG direction. In hgure 7 the blue solid (dashed) curve shows the expected anisotropy 
a-y (without) taking into account the absorption, for the hducial choice of lifetime discussed 
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Figure 7. The dipolar anisotropy induces by the 7 -rays, as defined in eq. (4.1), as function of energy. 
The blue solid (dashed) curve depict the anisotropy by taking into account (neglecting) the absorption 
of 7 -rays, for tdm = 10^® s. The red dot-dashed curve shows the anisotropy for tdm = 2.5 x 10^^ s 
which is the lower limit on lifetime at 2a from anisotropy data. The data points show the measured 
anisotropies by EAS-TOP [36-38], Akeno [39], IceTop [40] and IceCube [41] experiments. 


previously; while the red dot-dashed curve corresponds to the limiting value when exceeds 
the measured a at 2(t. For comparison, we also report the amplitudes of dipolar anisotropies 
measured by different experiments. A few remarks are in order: 

• The suppression of anisotropy due to the absorption can be clearly seen. It also con¬ 
tributes to the peculiar energy dependence of decreasing with energy, while the 
observed anisotropy a moderately increases with energy. 

• perhaps surprisingly, the bounds following from anisotropy are at least comparable 
in strength with the previously obtained bounds coming from comparisons with the 
(prompt) flux limits from EAS detectors and Fermi-LAT diffuse isotropic data, at the 
level of 10^^ s. 

• The observable, on the other hand, has a higher sensitivity to the inner Galaxy 
DM profile. For instance, the previously quoted bound of 2.5 x 10^^ s for the fiducial 
NFW profile would degrade to 1.9 x 10^® s for a cored isothermal profile [42] with 
Phio) = Ph/i^ + {p/f'cf') where Tc = 4.38kpc and = 1.387GeVcrn”^. 

One may wonder to what extent the above considerations are spoiled by a realistic 
account of experimental angular and energy resolutions. For instance, GASA-MIA had an 
angular resolution going from > 2° at low energies to better than 0.4° at high-energies [43]. 
This is almost irrelevant for a cored DM profile, while a 0(1°) resolution can degrade the 
bound for a NFW profile down to 6 x 10^® s. Nonetheless, this constitutes a sub-leading 
uncertainty with respect to the one coming from the unknown shape of the DM profile in 
the inner Galaxy. Goncerning energy resolution, it was demonstrated that GASA-MIA could 
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detect spectral features not larger than 0.2 — 0.3 dex in energy [44]. Recent cross-calibrations 
between the energy scale retrieved via surface detectors and the one inferred via fluorescence 
light suggest that there might be an over-estimate on the absolute energy scale of the former 
experiments of the order of 30%, see for instance [45]. Fortunately, despite these uncertainties, 
figure 7 shows that the expected signal peak is very broad, changing by no more than ~ 20% 
between 150 TeV and 350 TeV. We do not expect thus that accounting for energy resolution 
and energy scale uncertainty would degrade our conclusions by more than 0(10%). It is also 
interesting that at few hundreds TeV the expected anisotropy from DDM matching IceCube 
observations may be only ~ one order of magnitude below the measured overall dipolar 
anisotropy, while at higher energies (~ PeV) the suppressed anisotropy is smaller by a factor 
of few, and its ratio to the charged cosmic ray signal is significantly less favorable. This 
suggests a first potential strategy to improve the constraints by using the energy-dependence 
and phase information of the anisotropy: although no deterministic prediction of the expected 
anisotropy due to charged cosmic rays is possible, one could calibrate a model for the stochastic 
phase fluctuation on the high-energy bins (say, 700 to 2000 TeV) where the charged cosmic ray 
contribution to a is definitely expected to dominate, and use the low energy band (around 
200-300 TeV), where the fractional contribution of DM is expected to be maximal, to put 
a constraint on the amplitude of a sub-leading contribution to the total anisotropy due to 
a-j,. The latter is characterized by a fixed direction (constant phase) and a specific energy 
dependence, very different from the competing charged cosmic ray anisotropy. 

Despite some room for improvement, with past data one cannot really expect order-of- 
magnitude gains in sensitivity. However, a yet more powerful way to improve over the present 
analysis would be to rely on the gamma/hadron discrimination possibility attained by the 
present generation of EAS detectors. In general, cuts based on the morphology of the shower 
(sometimes dubbed “compactness” criteria) allow one to select a photon-rich sample, keeping 
fraction of the initial photons, while retaining only ^ £7 of the contaminating hadronic 
background. The ratio Q = allows to quantify the gain in sensitivity to a photon 

signal when this cut is applied. While some rejection capability was already present in old 
experiments, even gamma-ray astrophysics oriented EAS detectors of the past generation, 
such as MILAGRO, were limited by €h ~ 0.05 — 0.1 [46] with corresponding Q-factors never 
much larger than ~ 2. On the other hand, the situation is significantly different already at 
currently operating water Cherenkov EAS gamma observatories, such as HAWC. Such an 
experiment has similar energy resolution performances as the above-mentioned ones (about 
40% at > 10 TeV [46]), and even better angular resolution, about 0.1° at > 10 TeV [46]. 
But the major improvement is in the rejection capability of the hadron background: at high- 
energies, 99.9% of the background is routinely rejected, but stringent cuts with eh — 10“'^ 
and 67 ~ 25% above 10 TeV, i.e. Q ~ 30, have already been reported [47], with even better 
performances that could be attained [48]. With an effective area, A^s, approaching rsj 10 ^ m^ 
at high energy and a field of view of AD ~ 2 sr, HAWC is expected to reach a sensitivity 
below the level of the IceCube diffuse neutrino flux, thus providing a unique constraint on 
the electromagnetic counterpart of the neutrino signal [47]. A high-energy photon-enriched 
sample in T = 1 year of HAWC data would consist of about 

e.yTAQ Aeff (E-y > 10 TeV) ~ 10^ events , (4-2) 

against a number of background CR events e/i TAD Aeff 4>cr {Ecr > 10 TeV) ~ 3 x 10®. This 
would certainly ease the measurement of the gamma-ray spectrum, as already noted in the 
past [12]. However, it is perhaps even more remarkable that the expected anisotropies of 
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0 ( 10 %) in the gamma-ray sample correspond to variations in gamma-ray connts of ~ 10 ^, 
as opposed to anisotropies in the CR backgronnd which are expected to be of ~ 10^ events. 
Pnt otherwise, in the gamma enriched sample the anisotropy should be fully dominated by the 
gamma-ray contribution and, with snch statistics, HAWC may provide a crncial test of the 
DM hypothesis throngh anisotropy stndies, besides spectral ones. 

The sitnation may be even more favorable with fntnre detectors. HiSCORE [49] has two 
orders of magnitnde larger effective area than HAWC at high-energy, bnt 1 ^ Q ^ 2 and thns 
is not ideal for this kind of measnrement, althongh it may still be nsefnl for complementary 
stndies [12]. However LHAASO [50], thanks to its optimized hadron rejection capability, 
wonld provide the nltimate sensitivity for this type of analysis: According to [51], thanks to 
the KM2A array (for detecting hadronic indnced mnons) snrronnding the 10^ m^ Cherenkov 
detector, at > 80TeV LHAASO wonld reach e-,, ~ 1 and Ch < 10“'^, which ensnres that 
observations wonld be essentially CR backgronnd-free. 

In snmmary, anisotropy data offer a complementary tool to constrain DDM contribntion 
to the gamma-ray flnx at 0.1 — 2 PeV energy. A simple analysis shows that cnrrent constraints 
from the normalization of the anisotropy are competitive with the other methods, and we 
sketched two possible approaches to improve npon them: with a reanalysis of cnrrent data, 
the addition of phase information conld already allow one to achieve sensitivity to a snb- 
leading DM indnced anisotropy. The optimal energy window for DM signal to CR noise 
appears to be at ~ 200 — 300 TeV, while higher energies shonld be dominated by the CR 
component, and conld be nsed to calibrate the dominating backgronnd anisotropy, whose 
phase is flnctnating with energy. A major improvement relies however on the cnrrent (HAWC) 
and forthcoming (LHAASO) generation of EAS gamma-ray detectors: thanks to their greatly 
enhanced photon/hadron rejection capabilities, already in HAWC the anisotropy signal shonld 
be dominated by the gamma-ray one, allowing for a hrst stringent test of any “Galactic based” 
model for the origin of a signihcant fraction of IceCnbe events. In LHAASO, we expect the 
sample to be essentially backgronnd-free, allowing for the nltimate test (or detailed stndies, 
in case of positive detection) of these models. Needless to say, these experiments have a great 
potential for heavy dark matter constraints, which has only recently started to be stndied. 


5 Conclusions 

The IceCnbe discovery of a PeV flnx of astrophysical nentrinos has several implications for 
high energy astrophysics and astroparticle physics, as proven by the broad commnnity whose 
attention has already triggered. If the sonrces of (a fraction of) these events are Galactic, this 
discovery paves the way to “Galactic 7 -ray PeVatrons”. In this regard it is timely to investigate 
in more detail the pecnliarities of VHE 7 -ray propagation in onr Galaxy. In fact, the yet 
nndisclosed 0.1 — 10 PeV 7 -ray window in Galactic astrophysics wonld be affected by the pair 
prodnction absorption. In this article, we discnssed some effects that this phenomenon has 
onto expected signals in extensive air shower (EAS) detectors. We selected the benchmark case 
of a continnons emission from the Galactic halo in a decaying dark matter (DDM) scenario, 
althongh most of onr resnlts apply mutatis mutandis to other sonrce distribntions. Onr choice 
was also motivated by the observation that, while the potential of these EAS instrnments 
for DM detection has sometimes been considered in the past, see for instance [52, 53], their 
potential for serendipitons DM discoveries in more nnconventional scenarios has passed mostly 
nnnoticed. 
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In this article, we have calculated the expected 7 -ray flux from DDM, with values for 
the mass and lifetime motivated by the IceCube flux, by considering both the prompt and IC 
scattering components. While the angular dependence of the prompt 7 -ray flux is robust (with 
the detailed features of its energy spectrum depending however on poorly known branching 
ratios, and only computable within the theoretical uncertainties unavoidably affecting PeV 
scale physics), the IC scattering component can potentially exhibit unusual an^ii/ar features. 
Indeed, the IC scattering profile depends on the properties of the magnetized halo of our 
Galaxy, which is typically poorly known at large vertical distances from Galactic disk. By 
simple modeling of the Galactic magnetized halo, we have calculated the expected IC flux 
from directions near Galactic poles, arguing that in the most optimistic case the IC flux 
can be enhanced, becoming comparable to the prompt flux from this direction. However, a 
relatively large halo magnetic field at high latitudes will suppress the IC flux significantly. 

Typical EAS bounds on the 7 -ray fraction in cosmic ray flux have been derived under 
the hypothesis of a isotropic flux. We argued that this approximation is untenable in the 
energy range of interest, even in the limit of an isotropically emitted flux, due to the direction 
(and energy) dependent optical depth of the Galactic sky. We quantified this observation 
showing that the effect is so relevant that the current constraints from KASCADE or CASA- 
MIA experiments-which naively would appear to constrain some DDM scenarios~are still at 
least several times too weak. Ideally, an exposure only marginally better than the above 
mentioned ones, but at an observatory located in the Southern hemisphere, would be much 
more promising in this respect. We also argued that anisotropy may offer an independent 
handle to constrain DDM (as well as other similar scenarios): the expected 7 -ray flux induces 
an anisotropy in the overall cosmic ray flux only a few times smaller at ~ 0(100) TeV (and 
about one order of magnitude at ~ PeV) than the current measurements of the dipolar 
anisotropy routinely performed in EAS experiments. Turning the argument around, existing 
data are already sensitive to DM lifetimes of 0(10^'^) s, only one order of magnitude away 
from the value needed to fit IceCube events (~ 10^® s), showing the power of anisotropy 
analyses and motivating an attempt to improve over the current situation. 

Some progress is expected from the experimental point of view. For instance, the IceTop 
facility at the top of the IceCube detector can look at the Southern sky. Unfortunately, 
while located at the South Pole, the GC is not in standard analyses involving the IceTop 
array: the IceCube detector plays the role of penetrating muon detector for IceTop facility, 
which requires that the axis of air shower should pass through the volume of IceCube. This 
requirement leads to a cut on the zenith angle of shower ~ 30° for IC40 configuration [54], 
with the possibility of increasing to ~ 45° for the whole IC 86 configuration. At the South 
Pole, the GC is located at zenith angle ~ 61°. Although the expected sensitivity of IceTop, 
after 5 years of data taking, is close to the CASA-MIA limit at ~ PeV (see figure 14 at [54]), 
due to the closeness of the field of view to GC, IceTop can still moderately improve the limits. 
Both for past and forthcoming data, we argued that a dedicated analysis might be sensitive 
to sub-leading anisotropies expected from DDM, notably if the shape of the anisotropy and 
its energy dependence, whose expectations are relatively well-known, are imposed a priori in 
the analysis. 

Eventually, however, we have argued that greatly improved photon/hadron discrimina¬ 
tion capabilities are needed for a decisive jump in the sensitivity in both existing and forthcom¬ 
ing experiments: The recently inaugurated HAWC observatory [52] located in Mexico (with 
latitude A ~ 19° and zenith angle cut Om ~ 45°) thanks to quality factors Q = e^jy/eh ~ 30 
could provide first crucial tests of the DDM scenario, not only via spectral studies (see [12]) 
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but particularly when adding angular information, as discussed in our article. A future ex¬ 
periment such as LHAASO [50], benehting from the KM2A array, is expected to provide a 
gamma-enriched data set which is almost background (CR)-free, paving the way to exquisite 
constraints or, in case of detection, to detailed studies of the spectrum and morphology of 
the signal. 

Dehnitely, the opening of the PeV astrophysical window may offer new opportunities 
for interesting multi-messenger studies, probably shedding light on intriguing astroparticle 
questions. As already noticed in the past, and as we further argued here, searches in this 
new window signihcantly beneht from EAS experiments. Note that, while the low energy 
part of the neutrino flux observed by IceCube (recently extended down to ~ 10 TeV [55]) 
can naturally receive contribution from Galactic sources, perhaps easing their identihcation, 
pinpointing the origin of the high energy part of the flux is more challenging. In that respect, 
EAS experiments appear a unique and powerful probe. For such a task, as illustrated in 
this article, rather than considering the Galactic-scale horizon imposed by the hnite optical 
depth as a limitation, we should perhaps reconsider it as an original opportunity to exploit 
the specihc capabilities of EAS detectors. 
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A Details of the IC flux calculation 


In this appendix we provide the details of the IG flux calculation, which have been used to 
compute the curves in hgure 5. As we discussed, the IG flux is given by eq. (3.5) and here we 
discuss calculation of the ingredients dne/dEie and Piq. 

The steady state spectrum of at position x in our Galaxy (we will in fact assume 
azimuthal symmetry in cylindrical coordinates) can be calculated by 


dne 


{Ee,x) 


Ph{x) /■™ dm /2 


"IDMTDM b{Ee,x) 
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dEL 


{E'^)hisiE,,E'^,x)dEl 
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where b is the energy loss coefficient and /difj is the function taking into account the diffusion 
of e^. Let us elaborate on each of these ingredients: 

Energy loss function b{E(,, x): The lose energy by two main processes: synchrotron 
radiation in the magnetic held of Galaxy and energy loss due to IG scattering on the pho¬ 
ton bath (SL+IR and GMB). Since both the magnetic held and SL+IR held are position 
dependent, the energy loss function also is position dependent. So, 
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where ax is the Thomson cross section, T^ = Ae'j/me and 7 = Eg/me', and n(e,x) is the 
SL+IR+CMB photons differential number density at position x. The energy loss due to 
synchrotron radiation can be calculated by 


b (E 


(A.4) 


where B is the magnitude of the total magnetic field in our Galaxy, consisting of regular and 
turbulent components For the regular magnetic field we adopt the following profile [56] 
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where Rq = 8.3 kpc, tb = 10 kpc, zb = 2 kpc and Bq = 4.78 /rG. For the halo (possibly 
turbulent) magnetic field we assume a uniform constant strength magnetic field. Figure 8 
shows the energy loss function b, as function of Eg, at GG (black solid). Sun position (blue 
solid) and at vertical distance z = 5 kpc from Galactic plane (at GG) for three different as¬ 
sumptions for i?halo = 0,1 and 2 fiG, respectively by solid, dashed and dot-dashed red curves. 
As can be seen by increasing the value of i?haio from 0 to 2 fiG, the energy loss coefficient at 
z = 5 kpc increases by about one order of magnitude, which justify the suppression of IG flux 
(black and green dashed curves in figure 5) for larger i?halo- Obviously, the effect of halo field 
is smaller at lower energies, since the synchrotron emission is the main mechanism of energy 
loss at higher energies. 

Diffusion halo function Idiff(Ee, E^, x): The diffusion halo function can be calculated 
by solving the diffusion-loss equation of in the Galaxy. To avoid repetition we skip report¬ 
ing the details of calculation, which is done according to the prescription reported in [57]. 
However, it is worth mentioning that at high energies which we are interested in this paper 
-fdiff — 1) and so the results of of IG flux reported in figure 5 are only marginally dependent 
on the diffusion halo function. Put otherwise, the approximation described in the main text 
is actually excellent. 

The other ingredient in the calculation of IG flux is the IG power Pjc (see eq. (3.5)) 
which can be decomposed into the IG power from each component of the photon bath; i.e.. 
Pic = Bjc? where Pjc is the IG power from the photon bath n* with i = GMB and SL+IR. 
The PjQ can be written as 

Pfc(^7>^e,x) = ^^ l[ dq 

47 ^ 

(A. 6 ) 

where 7 = E(,/me-, e = E^jEf. and 
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^The interstellar magnetic field of our Galaxy off the Galactic plane is determined by the Faraday-rotation 
measurement of the polarization of extragalactic radio sources. However, these measurements depend on the 
assumed free-electron distribution at high latitudes which is poorly known. On the other hand, even in the 
Galactic plane, the estimated turbulent magnetic field is a factor of few larger than the regular magnetic field. 
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Figure 8. Plot of E~^b{Ee,r, z) as function of E^ at some positions in the Galaxy. For all curves, 
the assumed disk regular magnetic field is eq. (A. 5); while for solid curves i?haio = 0, for dashed red 
curve Bhaio = 1 mG and for dot-dashed red curve Bhaio = 2 /rG. 


Figures 9a and 9b, respectively, show the IC powers and as function of for 

Pg = 1) 10,10^, 10^, 10^ TeV at GC. As can be seen, at high energies the IC power sharply 
peaks at Eg. Namely, in the IC scattering almost all the energy transfers to the photon. 
Also, by comparing the corresponding curves in hgures 9a and 9b, it can be seen that at high 
energies the main contribution to the total IC power comes from The IC power in 

hgure 9a is independent of the position in Galaxy (due to the uniform CMB photon bath), 
while the IC power due to SL+IR strongly decreases by distancing from the GC, especially 
in the vertical direction with respect to Galactic plane. 




Figure 9. Inverse Compton power of: (a) P^^^{E^,Ee), (b) = 0), for Pg = 

1,10,10^, 10^, 10"* TeV, from the left to the right, respectively. 
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